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ABSTRACT 

A computerized technique combining backward wave propagation and an 
automatic edge detection scheme is developed and tested. The class of 
objects considered is limited to those with edge boundaries since it is 
shown that a universal automatic reconstruction cannot be obtained for 
all possible objects. Using samples of the acoustic diffraction pattern 
as input this technique enables the computer to predict the most likely 
locations of objects and to produce graphical output of the objects. 

A discussion is given on implementing the diffraction equations on the 
computer. A simplified edge detection scheme conserving both memory 
space and computer time is described. Test results are presented for 
both computer generated diffraction patterns and one set of experimental 
data. 
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I. MOTIVATION 



Since Dennis Gabor [Ref. 1] discovered holography - recording the 
amplitude and phase of the diffraction pattern of an object with the 
help of a reference wave and the subsequent reconstruction of the 
original object from this hologram - vast advances have been made in 
this field. Originally confined to optics holography has been extended 
to the fields of acoustics and microwaves; comparisons between the 
results obtained by the pioneers in holography and those of recent 
vintage are striking with the promise of even better results to come. 

(For a description of the development of holography see Gabor* s lecture 
on the occasion of receiving the 1971 Nobel Prize for Physics in 
Stockholm, Dec. 13, 1971, reprinted in [Ref. 2]). 

The interest in acoustical holography is motivated by the search 
for a means of probing media which cannot be investigated by any other 
method. For example, electro-magnetic waves have a very short propaga- 
tion distance under water (those that can be propagated for a longer 
distance at very low frequency, have much too low a resolution to be 
of any significance for imaging). In murky water optical wavelengths 
cannot be used at all. In another application it has been found that 
X-rays which give a clear indication of the bony structure of a body 
are not very good at making distinction between two different kinds of 
soft tissue (e.g. , in obstetrics). Acoustic waves however can show a 
difference in absorption for different soft tissues and do not result 
in damage to the insonified body at reasonable intensities. Investiga- 
tion by insonif ication and interpretation of the reflected sound has 
long been used in probing the earth* s crust and in nondestructive testing. 
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The application of holographic concepts brings to these fields the 
possibility of obtaining information that not only indicates the existence 
of an object of interest but also shows many of its features. The same 
statement holds for the advantages of acoustical holography versus the 
well known echo-ranging concepts as expressed in sonar. The result of 
the holographic approach is an image of the object, whereas without 
holography only detection of the object is realized. This is not meant 
to imply that the holographic approach is the only one feasible: 

Imaging systems incorporating acoustic lenses have also been developed 
with good results. (For the development and state of the art of 
acoustical holography the reader is referred to [Ref. 3]). 

In general illumination of a recorded hologram with a duplicate of 
the reference wave results in two reconstructed object wavefronts, one 
being an exact duplicate of the original object wave (producing the 
"true" image), while the other wavefront is the complex conjugate of 
the object wave (the n conjugate M image). If the two obtained wavefronts 
are not spatially separated from each other by some suitable arrangement, 
they can interfere with each other which results in distortion of the 
images (the so-called twin-image problem) . 

Unfortunately it is not possible to realize in acoustical holography 
that most striking feature of optical holograms, namely the three- 
dimensional sensation experienced by observing the true image with both 
eyes. In order to be able to retain the depth information, the wave- 
lengths of constructing and reconstructing beam have to be equal or 
nearly so. A hologram can be considered to be a sum of zone plates 
into which the object information is incorporated. (See [Ref. 4]). It 
then has an equivalent focal length f which is given by 
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f = 




-1 



where v and u are object and reference source distance respectively, 
measured from the hologram plane along the line joining object and 
reference source. See Fig. 1-1. 

Then if the linear dimension of the hologram is changed by a factor 
m from X to X 1 , (i.e., X 1= mX) , and the wave length X 1 used in recon- 
struction differs from the original wavelength X by y = X 1 /X , the 
effective equivalent focal length of the hologram is given by 



f 1 - ± U-^> 



-1 



where f}. = y*f/m 2 

v 1 and u 1 are the distance of reconstructed image and illuminating 
source respectively, measured as before, and the 4- sign gives the 
position of the true image, the - sign gives the position of the 
conjugate image. This gives rise to axial magnification 

M = (v/m 2 ) (l/v 2 ) 

A {+ (y/m 2 ) [(l/v) - (1/u)] + 1/u 1 } 2 

and transverse magnification 

These magnifications are different from one another unless m = 1/y 

and thus give rise to distortion of the reconstructed image, if this 

condition is not met (as is usually the case). In acoustical holography 

the reconstruction wavelength is optical. The sound wavelengths used 

in the recording process are much longer and y then is of the order of 
-3 

10 . To avoid distortion this would mean reducing the size of an 

acoustical hologram by this same factor, i.e., a hologram of lm 
diameter would be reduced to 1mm diameter. Naturally one cannot 
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expect to experience any three-dimensional illusion in looking at a 
hologram of this small size. Optical magnification of this small 
hologram introduces longitudinal distortion exactly analogous to that 
being eliminated by the scaling due to depth of field problems. This 
longitudinal distortion is one of the primary disadvantages of acousti- 
cal holography. 

Another disadvantage encountered in acoustical holography is the 
absence to date, inspite of intensive research, of a recording medium 
which gives the analogous fine grain structure as that encountered in 
optical holography by using spectroscopic plates. This shortcoming 
is in part balanced by the possibility of recording both amplitude 
and phase of a sound wave impinging on a linear detector directly. 
Because of the presence of linear acoustic detectors (one of the 
fundamental differences of acoustical holography versus optical 
holography) , it is not necessary to use an acoustic reference beam and 
to work with intensities as in optical holography. It might be argued 
that a direct recording of amplitude and phase of a diffraction pattern 
is not really a hologram (according to definition) because of the lack 
of the reference beam. . But as a hologram is after all an ingenious 
method of recording a complex wavefront using a medium that is insensi- 
tive to phase and as the recorded diffraction pattern contains the same 
information as a hologram, only in simpler form, there is no difference 
in the information content of the two types of recording. Working with 
amplitude and phase can also eliminate the twin-image problem referred 
to earlier. 

Assuming then the existence of a record of the amplitude and phase 

of the sound diffraction pattern due to some source (self-emitting or 
reflecting) , in sampled form, whether obtained by means of an array of 
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detectors or by scanning the sound field or by any other means, it is 
possible to circumvent the difficulties of optical reconstruction (not 
only depth distortion, but also those due to degradation because of 
nonlinearities of the film or cathode ray tube, etc. as in [Ref. 5]) by 
using a digital computer technique to restore the original object. 

However there will be errors introduced by quantization, sampling of 
the complex amplitude and nonuniformities in the scan if such is used. 
These effects are analyzed in [Ref. 5]. Of course the digital computer 
is limited in the amount of data it can handle, whereas the optical 
hologram is comparatively unlimited in the amount of information stored. 
Thus in increasing the number of independent samples of the wavefront 
eventually a point is reached in computer reconstruction beyond which 
the amount of storage and computation time becomes unacceptable. Another 
disadvantage of digital reconstruction is that perfect reconstruction 
can only be achieved for planar objects. The results of trying to 
reconstruct a three-dimensional acoustical object are identical to 
observing the real image in a reconstruction of an optical hologram. 

The shape of a three-dimensional object can only be reconstructed in a 
series of sections through its real image. 

One of the greatest drawbacks of digital reconstruction of objects 
from their diffraction pattern is that the distance from hologram to 
its generating object has to be known. This problem never occurs in 
the case of visual reconstruction of an object from its virtual image 
because the focus is found automatically by an eye - brain interaction. 

It is the same problem as that experienced in attempting to bring a 
real image to a focus on a screen, with the difference that the move- 
ment of the screen can be continuous and information on the validity 
of the placement is fed back by looking at the result and interpreting 
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it to be more or less correct than before. The final screen location 
is determined by the compatibility of the retinal image with some 
object hypothesis stored in the brain* s memory. In the field of eye- 
brain interaction great amounts of research effort have been and still 
are expended; still the processes are not completely understood to 
date. It is therefore hardly promising to try and simulate this inter- 
action on the digital computer. Nevertheless it would be of great help 
to find a method which brings about automatic focusing - or at least 
narrows down that region in space in which the object might be present. 
Problems can be envisioned in which acoustical holograms from quite a 
large number of different objects might be taken. Then one would not 
want to occupy a large number of specialists in the given field by 
having them look at the results of each backward propagation step but 
would rather have the computer make a preliminary decision on the object 
location and present this together with the computed sound field at 
that location. A single observer could then - using his pattern 
recognition ability - verify the computer decision or discard it. 
Naturally, computer time being expensive, the economics of either 
scheme have to be considered. In any case the computer program should 
take as little time and/or memory capability as possible, an obvious 
requirement for any efficient algorithm. 

This then is the problem to be solved in this investigation: Given 

a sampled sound diffraction pattern it was required to find a method 
for a digital computer to make a tentative decision on the object 
location. In all that follows a diffraction pattern will be assumed 



1A 



(except where stated otherwise) , that results from propagation through 
a medium that is linear and not disturbed by turbulence or convection. 
In the following section the necessary theory for this investigation 
is introduced. 



\ 



15 



II. SCALAR THEORY OF DIFFRACTION 



Diffraction is generally treated with a scalar wave theory. This 
is quite applicable for acoustical disturbances, more so than for 
electromagnetic waves which should rigorously be considered in a 
vectorial theory of diffraction as their electric and magnetic vector 
components are not independent but coupled through Maxwell* s equations. 

A. DIFFRACTION BY A PLANE SCREEN 

The solution of the problem of diffraction by a plane screen is 
treated in numerous books, e.g. ,. iRef. 6]. As it in itself is not a 
suject of this investigation only those results which are essential for 
the further treatment are stated below. Figure 2-1 shows the geometry 
of the problem. 

The resulting complex amplitude U (P Q ) at an observation point P Q 
behind an impenetrable screen S, with an aperture A, which is insonified 
by a disturbance JJ can be written as 

U (P Q ) = Op-G - U -p ) ds (2-1) 

4.7T “ 3 n 

S 1 

where (2 is the so-called Green's Function, which is chosen here to be 
a spherical wave expanding about the point Pq and ^/^n is the derivative 
in the direction of the outward normal to the plane of S. 

The difficulty encountered now is that the boundary conditions 
which should be applied to the solution of Eq. (2-1) are never known 
exactly, for even though the screen is assumed to be impenetrable to 
sound (except at its opening A), diffraction is known to take place, 
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Figure 2-1 

Geometry of Diffraction by a Plane Screen 
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the exact expression for which is to be found from Eq. (2-1). To 
simplify the problem, Kirchhoff set up the following conditions on S-^ : 

1. Across the aperture A the field distribution U and 3U/3h are 
the same as they would be without the screen. 

2. Over that portion of S-^ that lies in the geometrical shadow. 

U and au/8n are zero. 

It is noted that these boundary conditions are only approximations. 

Still Kirchhoff T s diffraction theory leads to very good results which 
are widely used in practice. 

Without going into the details (see [Ref. 6]) Kirchhoff* s diffraction 
equation for the plane screen insonified by a unity amplitude spherical 
wave emanating from 



where (n,rQ^) indicates the angle between the indicated vectors. 

Another severe inconsistency of the Kirchhoff boundary conditions is 
that by assuming a value of zero for both U and 3U/3n the problem is 
over-specified and requires the field to be zero over all space which 
clearly is not the case. Considerable criticism has therefore been 
raised against Kirchhoff **s treatment. Sommerfeld removed this diffi- 
culty by assuming a Green’s function which makes or 3G/3h identically 
zero over all of , without invalidating the development which leads 
to Eq. (2-1) (for.details see {Ref. 7J). His result is incorporated 
into the Rayleigh- Sommerfeld diffraction formula: 



( 2 - 2 ) 




exp [jk(r 2 i+ r 01 ) 3 

r 21 • r 01 

A 



[cos(fi,roi) - cos(n,Y 21 ) ]ds 



2 . < p o> ‘ 31 




(2-3) 



A 
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where U (P^) is the field distribution in the aperture plane. This can 
be expanded into (2-4) 



£ 



z o z i 






ff u (Pi) 

A 



exp(jk((x n -x 1 ) 2 +(y n -y 1 ) 2 +(Z f) -Z 1 ) 2 ) 1 / 2 )^ 
(XQ-Xi) 2 +(y 0 -yi) 2 +(Z 0 -Zj ) 2 



which shows that diffraction is a space-invariant phenomenon as the 
result depends on the difference between observation point and the 
point at which the disturbance is located. 

Diffraction also is a linear phenomenon, because the response to a 
sum of inputs is equal to the sum of the individual responses to each 
input, (within the limitations of this study). The linearity and 
invariance of diffraction imply that the result Uq due to some distur- 
bance U impinging on a diffracting aperture can be written in form of a 
convolution integral 

Uo = u * h 

where h is the so-called impulse response of the diffracting system. 

It is also possible then to write: 

F {1^} = F {U *h} = F {U} • F {h} = F {U}* H 



where F is the symbol for the Fourier transform, H is the transfer 
function of the diffracting system and use has been made of the con- 
volution theorem. 

For a unity amplitude expanding spherical wave incident on the 
screen Eq. (2-3) becomes 



U (P ) 



Jx 



ff 



exp[jk(r 21 +r 01 ) ] 



" 21‘ r 01 



cos(fi,r 01 ) ds 



(2-5) 



which differs from Kirchhoff's result only in the so-called obliquity 
factor cos(n,r Q1 ) compared to 1/2 [cos(h,T 01 ) - cos(n,T 21 )]. 
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The importance of Eq. (2-3) is that it allows the computation of the 
field at any observation point if the sound distribution in a plane 
across the direction of propagation is known. 

B. SPATIAL FREQUENCY APPROACH TO DIFFRACTION BY PLANE SCREEN 

The following discussion follows closely [Ref. 7 and Ref. 8]. 
Propagation and diffraction being linear, space invariant phenomena it 
is possible to define impulse responses and transfer functions for 
different systems. Every wave distribution can be expanded into a 
finite or infinite sum of weighted plane waves. This is accomplished 
by Fourier transforming the complex disturbance by use of the formula 

00 

A (f x ,f y ) = J7U (x,y) exp[-j2TT(f x *x+ f y *y)] dx dy (2-6) 

— 00 

where A (f x ,fy) is then the complex spatial frequency spectrum of IJ 
and the spatial frequencies f x and fy are related to the direction 
cosines by 

f x = ct/A ; f y = 6/A; f 2 = y/A 

The direction cosines a, 8 and y are defined to be the cosines of the 
angles between the direction of propagation of the plane wave and 
the positive coordinate axes. The equation for a unit-amplitude plane 
wave propagating with direction cosines (a,8,y) is simply 

B. (x,y , z) = exp [jk(ax +8 y +y z)] (2-7) 

where the direction cosines are geometrically constrained to obey the 
equation a 2 + 8 2 + y 2 = 1 

Writing IJ (x,y,z) as an inverse transform of its spectrum 
00 

U (x,y,z) = ff A (f x ,f y ) exp[j27T(f x «x + f y *y)) df x df y (2-8) 
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it is obvious upon comparing Eq. (2-8) and Eq. (2-7) that U_ can be 
visualized as a sum of plane waves propagating with direction cosines 
a = X-f x , 6 = X *f y y = (l-(X 2 f x 2 )-(X 2 f y 2 )) 1/2 

The complex amplitude of a plane wave component is simply A (fx,fy) 
df x dfy, evaluated at (f x - oi/X,fy = (3/X) . 

In the following treatment the impulse response and transfer function 
of a diffracting aperture will be found. An equivalent form [Ref. 8] of 
the Rayleigh - Sommerfeld diffraction formula Eq. (2-3) is 



I (».) - h ff 3L CP,) ds 

A 1 r 01 



(2-9) 



= 27 // "<*!>• *01 ds 

A 

3 exp(jkr 01 ) 
where Ifo- ai” [ ] 

The plane wave expansion for a diverging spherical wave of the form 
[exp(jkrQ^) ]/ tq^ is known to be (see [Ref. 9]) 



exp(jkr Q1 ) j 1 

=— //y exp(jk[(x 0 -x 1 )o+(y 0 -y 1 )e 



01 



(2-10) 



+ ( 2 q- z i)y]> dadg 

To find the differentiation according to z^is performed and yields: 



2tt. 



Kqi = — ' ff exp (jk[ (x 0 -xi)o+(y 0 -y 1 )B+(zo-z 1 )Y] } dctdB (2-11) 



Then 



jj (V 



— ff u. (Pi) •// exp{ jk[ (x 0 -x 1 )a+(y 0 -y 1 )6 
X A 

+ [ (zq-z^y ] }dadB dx 1 dy 1 



( 2 - 12 ) 
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(2-13) 



This can be brought^into the form 

— (x 0 ,y 0 » z o ) = ^ X 1 ’ ^1 ’ Z P ' ^ ^ X 1 ’^1 ’ Z P 



a 3 



• ff exp{j2|{(x 0 -x 1 ) a +(yo-yi)3-t-(zo' z lH]} d X d Y 



by introducing the aperture function t(x^,y^ 5 z^) which is unity inside 
the aperture located in the opaque screen at and zero outside A. 

Kirchhoff boundary conditions are applied hereto U but not to 3U/9n. 
Separating the exponential and interchanging the order of integration 
leads to 



ot $ 

Uo(x 0 ,y 0 , z o) " ff ff UtCx^yi.z-p expHauO^-tyy-^] dxjdy-L 



— CO —00 



•exp{-j-^-yzi} exp{ j|^(otx 0 +3y 0 +Yz 0 )}. d^ d| 



(2-14) 



where Uj. in Eq. (2-14) replaces Uj_*t. But now the inner integral is 
just the Fourier transform or plane wave decomposition A^of U^ 



a 6 



At = ff UtU-py-^zp exp[-j2TT(yx 1 +^-y 1 )] dx-jdy-j^ 



(2-15) 



and the outer integral can be written as 

°° B TT 

Uo(xo,yo,z 0 ) = ff APX’X^ ex P^X“ (Z 0 “ Zp} 



* exp{j2^.(ax 0 +3yo)><^- d~- 



(2-16) 



a 6 

Eq. (2-16) can be seen as an inverse transform which relates Aovrjf*) 

^ A A 

to ^(xQ.yQ.Zo^ ’ so that 

4^. (— >— ) exp [ jZl-(zg-z^) /I-aZ-3^1 . must be the spectrum AnO^-rf) 

X X x A A 

of the diffracted and propagated disturbance. Therefore a transfer 

function H (f ,f ) can be defined 
— v x* y' 



« - ^Slcv ) ‘ oxp[j2l(z 0 - Zl )/l-(X£ x )^(Xf y ^] (2-17) 
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The corresponding impulse response is found directly from the Rayleigh- 



Sommerfeld integral Eq . (2-3). 

1 exp(jkr 01 ) ^ _ 

U < P o> = j X//1 ( p i) ^ cos (n, r 01 ) ds = 

A 

• i 00 exp(jkr n1 ) A - 

= — If U (P-. ) *t(Pi) cos(n,r 01 ) ds 

jA — L 1 r 01. 

oo 

CO 

= ff Ut( p i)‘h C p 0» p i) ds 

CO 



The impulse response Ti (Pq.,P^) is given thus by 

h (Po.P!) - 4 M * Uk, M 1 » .(a.7 01 ) 

“ jA r 01 01 

By definition then Eq. (2-17) must be the Fourier transform of 
Eq. (2-18). 



(2-18) 



(2-19) 
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III. FRESNEL DIFFRACTION AND DIFFRACTION 

TREATMENT IN THE SPATIAL FREQUENCY DOMAIN 

Having found the general form of the solution to diffraction by 
a plane screen it is soon evident that the Rayleigh-Sommerf eld 
diffraction formula is not very amenable to closed form solution. 

Even on the digital computer one would expect quite long computation 
times. Therefore approximations are introduced, which will allow 
reducing diffraction-pattern calculations to simpler mathematical 
operations . 

A. FRESNEL DIFFRACTION 

Consider the situation of Fig. 3-1. Single-frequency sound is 
diffracted by a finite aperture A in an infinite plane opaque screen 
which is located in the x^-y^- plane at the origin of the z- axis. 

The observation region is assumed to be parallel to the x^-y^- plane 
and has a coordinate system (xg , yg ) attached to it. This observation 
plane is at a distance Z from the diffracting screen. 

Using the Rayleigh-Sommerf eld formula the field at a point 
(x ,y ) can immediately be written as 

00 

Uo( x 0 > y o» z o^ = H - (VWV ^t* x i ,y i*°* dx i dy i (3-D 

—00 

where 

1 exp(jkr 0 i) 

h. ( x 0» y 0‘ x l ,y l^ = TT r cos (n,r Q1 ) (3-la) 

01 

and the effect of the finite aperture A is incorporated into U it 
being assumed the U^. is zero outside the aperture region. This allows 
extending the integration limits to infinity. 
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The first assumption to be made is that the distance z is very much 
greater than the maximum linear dimension of the aperture A and that in 
the observation plane only a small region around the z-axis is of 
interest. Again z is assumed to be very much greater than the extension 
of this region of interest. The initial approximations introduced are 
called the paraxial approximations because they hold only close to the 
axis of propagation. Using this assumption the obliquity factor is 
approximated by 



which is good to within 5% if the argument of the cosine is not greater 
than 18°. 

The same assumption can also be used to replace the radius by 

the distance z in the denominator of Eq. (3-la) but not in the exponent, 
since the resulting error would be multiplied by the large number k, 
generating intolerable phase errors. The impulse response has now been 
approximated by: 



Further simplification is possible by using the Fresnel approxima- 
tion. This consists of expanding the radius r^ in the exponent, which 
is given exactly by 



cos (n,r 01 ) Hi 



h ( x 0 »y 0 » x 1 »y 1 ) ~ exp(jkr oi ) 



J AZ 



(3-2) 



r, 



01 = Ax 0 -Xi) z + ( y 0 -yi)‘ i + zZ 



z /l+( Xo-xi^ + ^yp-y fy 2 



Z 



Z 



in a binomial expansion and assuming that r^ is given adequately by 
the first two terras in the series. Then 



r 



01 s Z[1 + 1 <i£til)2 + A («£&)* J 
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The impulse response can then be written as 



h (x 0> y 0 ,x 1 ,y 1 ) = e - X -^ 2" Z - ) - ex P ^2? [(x 0 _x l )2+(y 0" y l )2;i} (3-3) 

The region for which this approximation is exact enough is called the 
region of Fresnel diffraction. The accuracy of this approximation 
hinges on the relative sizes of the aperture and of the observation 
region, the propagation distance z and the wavelength X. One might 
require that the maximum phase change contributed by the third order 
term in the series expansion should be rather small, say much less 
than one radian. This is the case if 

Z3>> 4l [(x 0“ x l )2+(y 0- y l )2] HAX (3-4)’ 

For an observation region and diffraction region of about 50X by 50X 
this would require the distance z to be much greater than about 250-300X. 
In fact the requirement of Eq. (3-4) is too restrictive because all 
that is required is that the higher order terms do not change the value 
of the superposition integral of Eq. (3-1) by an appreciable amount. 

For distances z smaller than those given by Eq. (3-4) the oscillations 
of the quadratic phase factor given by Eq. (3-3) will generally be so 
very rapid as the integration is performed that they virtually cancel 
out and only those points of the aperture where Xq=x^ and ^0 = ^1 
give a noticeable contribution to the integral. Near these points of 
so-called stationary phase the magnitude of higher order terms will 
normally be negligible. The principle of stationary phase is discussed 
in [Ref. 6]. None the less, Fresnel diffraction cannot be assumed to 
hold for distances which for the given example might be less than about 
100X (the paraxial approximation is implied also !) . 
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If one is satisfied as to the validity of the Fresnel approximation 
for a given case, the superposition integral can be expressed in' the 
form of the following equation 



exp(jkz) 



J k 



U (x 0 ,y 0 ) = ■ exp ^ (x 0 2+y 0 2)] 



J A 2 



(3-5) 



r /„ 2 i 2> ’ - .2ir 



* si u t (x i ,y i^ exp ^Tz ( x i + yi ex P I-jxz ^ x o x i +y o y i^ 

; . dx x d yi 



by expanding the quadratic terms in the exponent. This can be seen 
to be (aside from multiplicative phase factors and attenuation indepen- 
dent of (x-^,y^)) the Fourier transform of U t (x^,y^)* exp[iii(x^ 2 +y^ 2 ) ] , 
evaluated at the spatial frequencies f^ = Xq/x z and f^ = yg/xz, which 
assures correct scaling in the observation plane. As such the integral 
Eq. (3-5) is amenable to digital computer solution using a Fast Fourier 
Transform algorithm. 

It is also possible to transform the impulse response of Eq. (3-3) 
because it is space-invariant (it depends only on the differences 
(xq-x-l) and (yg-yi)) to get a transfer function: 



H (f x »f y ) = exp(jkz) exp[-jTrAz(f x 2 +f y 2 )] , (3-6) 

which relates the spatial frequencies in the diffracting and in the 
observation plane. This is immediately observed to be an approximation 
to the exact transfer function of Eq. (2-17). 
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B. FRAUNHOFER DIFFRACTION 



For more severe restrictions on the propagation distance, it is 
possible to simplify the diffraction equation even more. If z is large 
enough to assume that the quadratic phase factor exp[ jk(x^2+y^2) /2z ] is 
approximately unity, then the diffraction equation reduces to a straight- 
forward Fourier transform (again evaluated at the scaled values of 
spatial frequency) with an additional multiplication 

U (x Q ,y 0 ) = ex P i J kz . ) - exp [jk(x 0 2+y 0 2)/2z] 

J A z 

“ (3-7) 

X U £ t ( x i»yi) ex P I-j 2 ^(x 0 x 1 +y 0 y 1 )/xz] dxjdyj^ 

—CO 

This equation is called the Fraunhofer diffraction formula. Unfortu- 
nately, the distance z would have to be very large indeed for Eq. (3-7) 
to be a good approximation to the actual physical situation. For the 
dimensions used previously to find the minimum distance for Fresnel 
diffraction z would have to be very much larger than 15000A. Neverthe- 
less, the Fraunhofer diffraction equation can be very useful in optics, 
because it fits the situation at the focal plane of a well-corrected 
lens; it is of less use. in acoustical diffraction problems because of 
the difficulty of building good acoustical lenses. In the context of 
this investigation it can be used to check the results obtained by use 
of the Fresnel equation because Fraunhofer diffraction is the limit 
of very large propagation distance z in the Fresnel diffraction equation 
and therefore a Fourier transform of the diffracting aperture would be 
ob tained. 
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The transfer function QSq.- (3-6)) approach to Fresnel diffraction can 
also be used in the Fraunhofer diffraction region, but one would hardly 
under normal computational circumstances try to implement the Fresnel 
(or Fraunhofer) diffraction equation by using the transfer function in 
the spatial frequency domain. This will take longer computation time 
by requiring a forward and an inverse Fourier transform, plus a matrix 
multiplication whereas in programming the Fresnel equation in the space 
domain one needs only a forward Fourier transform with additional 
matrix multiplications, 

C, SPATIAL FREQUENCY APPROACH 

If it has been decided that the use of longer computation time is 
feasible, then the implementation of the spatial frequency diffraction 
equation (Eq. (2-26)) is a much more powerful tool. (Indeed, for 
large array sizes the use of the spatial frequency approach may be more 
economical, as the number of computations for a Fast Fourier Transform 
increases like N*log 2 N, whereas N 2 multiplications are required for each 
of the quadratic phase factors). By using it not only the Fresnel 
requirement on propagation distance but also the paraxial approximations 
are circumvented. Any inaccuracies in the results obtained by the 
spatial frequency approach are thus due only to the use of Kirchhoff 
boundary conditions on JJ or on 3U/an. But as Kirchhoff T s diffraction 
theory as used in the Fresnel-Kirchhof f diffraction formula or in the 
Rayleigh-Sommerfeld diffraction formulation (the relative accuracy of 
both equations being still a subject of active research) gives results 
which are in excellent agreement with practice, any reservation against 
the use of Kirchhoff’s boundary conditions can only be very slight. 
Moreover they are applied in the Fresnel equation also. 
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To program diffraction by a plane screen using the spatial frequency 
approach one first multiplies the insonifying disturbance by the aperture 
function, performs a forward Fourier transform and multiplies each 
spatial frequency component by the appropriate phase delay as given by 
Eq. (2-26). Inverse transforming then results in the diffraction 
pattern in the observation plane. 

In shorthand notation: 

U (x 0 ,y 0 ) « F-l {H (f x ,f y ) • F {UjlCx!, yi ) • t (x 1>yi )» (3-8) 

where F and symbolize forward and inverse Fourier transform and 
— ^ f x ,f y) is § iven by Eq ‘ ( 2 “ 2b ) 

H (f x , f y ) = exp tj^p- A-CXf^^-Ufy)^] (2-26) 

1J^ is the insonifying wave complex amplitude and t_ is the aperture 
function. 

The aperture function is usually taken to be unity inside the 
aperture and zero outside but there is really no necessity for doing 
that. The aperture might also be characterized by the introduction e.g., 
of a constant phase shift or a quadratic phase shift (like that due to 
a lens) or a combination of amplitude attenuation and phase shifts. 

Of course the effect of the aperture could also be handled in the 
spatial frequency domain by using the convolution theorem of linear 
systems analysis and writing 



F •t.( x 1 >y 1 ) ) ■ A^x.fy) * 1 ( f x > f y ) 

where A (f x ,f^) is the spatial frequency spectrum of UpT^(f x ,fy) is 

the spatial frequency spectrum of _t and * represents the convolution 
operation. 



31 



This explanation may be helpful in visualizing the effect of the 
aperture in broadening the angular spectrum of the disturbance but would 
not be used on the computer if it can be avoided. 

The transfer function H (f , f ) has to be investigated more closely. 

x y 

To recapitulate 

H (f x ,f y ) = exp A-Uf x )Mxf y )*] 



(2-26) 



When the direction cosines a=X*f and 3 =X • f satisfy a 2 +3 2 <l (case of 

x y 

homogenous waves) the effect of propagation over a distance z is just 
a change in the relative phases of the spectral components: Since each 

component travels at a different direction it has to travel a different 
distance to reach the observation plane and thereby relative phase shifts 
are introduced. But when a 2 +3 2 >1 (case of inhoraogenous or evanescent 
waves) the square root in Eq, (2-26) is imaginary and 

H (f x ,f y ) = exp Af x ) z + (xfyT^l] = exp I-pz] 

where p is a positive real number. These wave components which travel 
perpendicular to the direction of propagation are severely attenuated 
by even small propagation distances. For propagation distances greater 
than a few wavelengths the influence of these evanescent waves is 
entirely negligible on the diffraction pattern and it is possible to 
band-limit the transfer function to 



H (f X ,£y) - 



exp[j-?15. /l-(Xf x ) 2 -(Xf )‘ ] for f 2 +f 2 < 4— 
X y x y x 2 

0 elsewhere 



(3-9) 
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D. COMPARISON OF FRESNEL EQUATION AND SPATIAL FREQUENCY APPROACH 

A comparison between the two possible ways of implementing the 
diffraction problem is now possible. 

The time needed for one propagation step using the Fresnel equation 
(one Fourier transform, two multiplication operations on a two-dimen- 
sional matrix - one before, one after the Fourier transform) will be 
less (for moderate array size) than that consumed by one propagation 
step using the spatial frequency approach (forward Fourier transform, 
matrix multiplication, inverse transform). The Fresnel equation 
cannot be used for small propagation distances because of its inherent 
approximations, which do not occur in the spatial frequency method. 

The spatial increment between array points increases with propagation 
distance in the use of the Fresnel equation due to the necessary 
evaluation of the Fourier transform at fx = x Q /^z, fy=y 0 /x z > w ^ ereas 

it is constant .when the spatial frequency approach is used. This leads 
to the fact that it is possible to use the Fresnel equation program for 
very large propagation distances whereas the spatial frequency method 
is limited in the propagation distances to be used because the diffrac- 
tion patterns from "ghost" objects (which are postulated by the Discrete 
Fourier Transform) start overlapping with the useful diffraction pattern. 
These effects are described in more detail later in the section on the 
Discrete Fourier Transform. Reduction of this unwanted and damaging 
effect is possible but only by using a larger array size which can 
increase the computation time considerably. The last aspect under 
which the two methods will be compared is that of the validity of the 
results after backward propagation of the observed wavefront when the 
propagation distance differs by a small amount from the actual distance. 
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Suppose a diffraction pattern is observed at a distance z from a plane 

object. Then, on propagating the wave back a- distance z-dz (an error 

which might occur in an automatic reconstruction scheme) using the 

transfer function approach one would expect to obtain the same 

diffraction pattern as would appear on forward propagation of the 

original diffracted wavefront at the object by a distance dz. It is 

then possible (by propagation backward using a propagation distance dz) 

to obtain the actual diffracting object. The case is quite different 

when the Fresnel equation is used. Here when the propagation distance 

in reconstructing differs by a small amount from the true one, the 

obtained "image" is only a very coarse approximation to the diffraction 

pattern at a distance dz in front of the diffracting aperture. The 

actual diffracting plane and that just reconstructed are related only 

by the fact that they give rise to the same diffraction pattern at a 

distance z and z-dz respectively when the Fresnel equation is used. 

In this second case it is not possible by another small propagation 

step of the correct difference dz to reach the correct result: The 

Fresnel equation does not hold true for small propagation distances. 

One would have to go back to the original diffraction pattern and to 

try again from there. This either requires storing the original 

diffraction pattern (because in the use of the Fourier transform the 

original data are replaced by the transformed components), which being 

a two-dimensional matrix of complex quantities with a rather large 

size requires a considerable amount of additional core memory or 

demands another propagation step for the restoration of the original 

pattern. Thus the use of the Fresnel equation in an automatic 

reconstruction scheme would take more computation time or core memory 

than the spatial frequency method and become uneconomical. (If the 
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propagation distance is known precisely beforehand this reversal 
of relative time consumption between the two methods would not occur) . 

This last aspect settled the question which algorithm should be used in 
favor of the transfer function approach in that it is desirable in a 
restoration problem to reach a result which is as close to the actual 
disturbance at a particular distance as possible even though the 
distance might be slightly wrong. Moreover it is much more difficult 
to handle a diffraction problem by the Fresnel equation if two diffracting 
objects are located in planes which are displaced in depth, due to the 
self-scaling of the Fresnel equation; it is impossible to do if the 
distance between these objects is less than that required for the 
Fresnel equation to hold. 

E. BACKWARD PROPAGATION 

How .'backward propagation' is accomplished can now be described. 

The conceptually more satisfying spatial frequency approach is treated 
first because it is used in the actual program. Consider a given 
diffraction pattern Uq (xQ,yg) at a distance z from an object H]_( x i>yi) 

as in Fig. 3-1. Then from the diffraction equations, 

Ho (*o.y<>) = F " 1 ( f x’V‘ F fui< x i»yi)» (3-8) 

To find ( x i»yp one would first perform a Fourier transform on Ho> 
multiply the result by H - ^ and then inverse transform the result 

Ui (x^yp = F" 1 (Uq ( x q>Yq) } } (3-10) 

where H“^ is defined by 

H-H _1 - 1 (3-11) 
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The symmetry between Eq. (3-8) and (3-10) is obvious. The reverse 
transfer function H” 1 can easily be found from the definition Eq. (3-11) 
and Eq. (3-9) which defines the transfer function H. 




(3-12) 



Observe that the same result would have been obtained by mechanically 

inserting (-z) as the propagation distance in the equation for the 

forward transfer function, where the minus-sign might be taken to 

signify backward propagation*. Another observation is in order here: 

Suppose II were not taken to be band-limited to those spatial frequencies 

which give a real square root. Then at least HT^ must be band-limited; 

\ 

else for those frequencies which give an imaginary square root 

if 1 (f x ,f y ) = exp [-^p Axf x ) 2 + (Xf y ) 2 -1] =e P for f x 2 +f y 2 >-^y 

and an exponential growth would be postulated, which leads to conceptual 
problems (for complete treatment of this subject see [Ref. 8]). Apart 
from these conceptual problems: the computer will not handle exponentials 
with real positive arguments which are too large (as would be in the case 
for non-bandlimi ted II ^ used in a propagation over more than a minuscule 
distance) . The frequency components of the angular spectrum correspond- 
ing to evanescent waves in a diffraction field that is several wave- 
lengths away from the object would be very small in amplitude. There 
is a good chance that they would be prone to round-off errors on the 
computer and in an actual data acquisition system would be lost in the 
unavoidable noise anyway. # Trying to work with a non-band limited 
reverse transfer function then is bound to fail. But it must always 
be kept in mind that the evanescent waves are physically existent and 
carry information about those object features which are smaller than 



36 



a wavelength. Setting them to zero destroys information irretrievably 
and therefore degrades system performance. Figures 3-2 through 3-5 
give an example of this image degradation. Figure 3-2 shows the 
amplitude distribution of the object, which is 32 sampling points wide 
and 16 sampling points deep, centered in a 64 sampling points square 
array. If the emission at all sampling points is assumed to be in 
phase (as in a specular object) and the sampling points are one half 
wavelength apart the result after transforming and reconstruction 
neglecting evanescent waves is given in Fig. 3-3. The degradation is 
not very severe. For in-phase emission and sampling every 0.3 wave- 
lengths the result of Fig. 3-4 is obtained which shows much lower slope 
and higher overshoots. It still gives a good indication of the actual 
shape of the original object. Figure 3-5 is the result of sampling 0.5 
wavelengths apart and assuming a random phase distribution for the 
object to be reconstructed (as in a diffuse object). As a diffusely 
emitting object will have more energy in high frequency components 
one would expect the loss of information to be worse than in the 
specular case; still the degree of degradation as shown in Fig. 3-5 is 
rather unexpected. For experimental results of backward propagation 
using the spatial frequency approach see [Ref. 10]. 

From the Fresnel equation U-^ (x^,y^) can be found by forming the 
complex conjugate of (x 0 ,y 0 ) propagating Uq* a distance (+z) 

,M x 2»y2) = - e tf Z) ex P ^2z (x 2 2+y 2 2) J (3-13) 

oo 

• // Uq* (x 0 ,y 0 ) exp{jj^(x 0 2 +y 0 ? ) ] exp [-j2l(x 0 x 2 +y 0 y 2 ) ] 

oo 2 Z Az 

dx o d y 0 
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Figure 3-2 

Amplitude Distribution of Input Object 
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Figure 3-3 

Reconstructed Specular Object with Nyquist Rate Sampling and Neglecting Evanescent Waves 
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Figure 3-4 

Reconstructed Specular Object with Sampling at Higher than Nyquist Rate and Neglecting Evanescent waves 
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Figure 3-5 

Reconstructed Diffuse Object with Nyquist Rate Sampling and Neglecting Evanescent Waves 



Conceptually, the term ’backward propagation’ is not applicable to this 
operation; it is rather a forward propagation of a different though 
closely related diffraction pattern. The result of this operation is 
the original object with the positive sense of the x^ and y^- axes 

reversed, i.e. , U2(x2,y2) - H]/"~ x l , ~ y l^ this usually is no reason 

for not using this approach, since the inversion is easily corrected. 
The main reason for discarding the Fresnel equation approach to back- 
ward propagation is that it is not possible to use repeated small 
propagation steps in the Fresnel equation due to its inherent 
approximations . 
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IV. DISCRETE AND FAST FOURIER TRANSFORM 



Diffraction pattern computations can be accomplished in various 
ways, but if one of the approaches as outlined in the last section is 
chosen (for the sake of simplicity and short computation time) it will 
be necessary to take the Fourier transform of some input array. Since 
it is possible to work only with discrete samples on the digital 
computer a discrete Fourier transform has to be used. 

A. DISCRETE FOURIER TRANSFORM 

Suppose a function f(x) is represented on the interval [0,X] by N 
samples f(n*a) where 




0 < n <_ N - 1 

The Discrete Fourier Transform (DFT) of the sequence f(o), f(a),*"* 
f((N-l)*a) is defined to be the sequence F(o), F(R),**‘ F((N-1)«0) 
where 

N-l 

F(kfl) = I f(n*a) exp(-jn'a*k«0) 0 k _< N-l (4-1) 

n=o 

and Q = — = h— (4-la) 

X N-a 

The inverse DFT is given by 
1 N-l 

f(n*a) = — I F(k* 0) exp(jkOna) (4-2) 

w k=o 

Remarks : 

1. The similarity between the DFT as defined above and the 
exponential form of the Fourier series expansion is obvious. There 
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f(t) = i E F(noi) exp(jntot) 

and t/2 

F(nm) = f(t) exp(-jnmt) dt 

-T/2 

In the case of the DFT the integral is replaced by a sum as there are 
only discrete samples of f(x). 

2. The DFT gives as many frequency components as there are samples 
of the function to be transformed. 

3. In order to make the definition of the DFT in form of a Fourier 
series consistent it has to be assumed that the sampled function is 
periodic with period X. Otherwise it could not give rise to discrete 
frequency components. Thus each object must be taken to have 'ghost 
objects' of equal shape located next to it at a periodic distance X 

in both horizontal and vertical directions. Analogously, the sequence 
of spectral components must be assumed to be periodic outside the 
actually computed spectral period. 

4. Even if the sampled function values are real, the sequence 
of frequency components is in general going to be complex. 

5. The spacing between frequency components is the reciprocal 

of the period X, multiplied by 2 tt. As the spatial frequency approach 
to diffraction was shown to be bandlimiting in the last section, it 
is possible to compute the minimum sampling rate, which is necessary 
to avoid overlapping of spectral components. In the spatial frequency 
domain the spectrum of the transformed object is repeated with a 
frequency period of 1/a. If a is chosen to be large the individual 
spectra will begin overlapping; it is then not possible to extract the 
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original object from this distorted spectrum. This effect is called 

'aliasing 1 . To recapitulate: 

e*p(j^P /l-(>f x ) 2 -(Xf ) 2 ) (f x 2+E v 2< “ 

H - * y X y A 2 

0 elsewhere 

This means that in the spatial frequency plane H (f ,f ) is 

x y 



bandlimited to those frequencies which lie inside the circle of radius 
l/X centered at the origin. (See Fig. 4-1) If efficient sampling 
is desired (meaning sampling that is not higher than the Nyquist rate) 
the sampling scheme becomes very complicated even for such a simple 
band-limited frequency region as a circle. This can be avoided at the 
cost of slightly higher sampling rate (about 12%, see [Ref. 11]) by 
assuming the band-limited region to be the smallest square that circum- 
scribes the actual circle. Using the Nyquist sampling criterion on 
this new band-limited region leads to the maximum allowable spacing 

t 

between sampling points, namely a =X/2. This is found the following 
way : 



For the f x -direction: 



f :xMAX X 



Af x N/2 




J2 

2tt 



from Eq. (4- la) q = — 

X 



then 



V _ NvX. 
X 2 ’ 




2-2tt 

N-X 



X. 

2 



A similar result holds analogously for the fy -direction. 
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Figure 4-1 

Propagation Transfer Function in the Spatial Frequency 
Domain, Showing Band-limitation at Nyquist Rate Sampling, 




Figure 4-2 

Propagation Transfer Function in the Spatial Frequency Domain, 
Showing Band-limitation at Twice Nyquist Rate Sampling. 
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Now, suppose for some reason it would be desirable to introduce a 
guard band in the spatial frequency domain or to sample at a higher than 



the Nyquist rate, e.g.. 



1 A . X 
a ’ 4 <a ~:1 

without increasing the total number of samples. 



Then f 



« 1.-2 



xHAX 2a 1 



= _ . (See Fig. 4-2) 



But the Discrete Fourier Transform will space the calculated frequency 
components evenly from -f x max to 4- f^ax, and from -f^max to +f y max 

and all the spectral components outside the bandlimiting circle will be 
set to zero on propagation. Therefore the (necessarily discrete) 
propagation algorithm will have to work on fewer frequency components 
than when the Nyquist rate is used. The results will be in error 
because on inverse transforming fewer spatial frequency components 
are available. There is then a trade-off between fine sampling in 
space and a coarsening of results of the propagation algorithm which 
has to be taken into account. 

Of course if the object spectrum itself is band-limited to a band 

which is narrower than that introduced by the propagation transfer 

function, the sampling rate could be adjusted accordingly to a coarser 

sampling (sec [Ref. 12]). Since the object is assumed to be unknown, 

this possibility cannot be used here. 

The definition of the Discrete Fourier Transform is readily extended 

to two (or more) dimensions. In particular for two dimensions: 

N-l M=1 

F (kft,lA) = £ E f(na,mb) exp [-j (nakft 4- mblA)] 

n=o m=o 

for 0 < k < N-l and 0 j< 1 <_ M-l (4-3) 

where ft - and A = a= — ; b “ ~ 

a Y N M 
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Then 



] N-l M-l 

f(na,mb) = Tfn7 £ l F(kft,lA) exp [ j (nakft+mblA) ] (4-4) 

" k=0 1=0 

The postulated periodicity of the object to be transformed is very 
unfortunate for a system that is to simulate imaging and diffraction. 
Consider the one- dimensional case as shown in Fig. 4-3. An object and 
its two nearest ’ghost objects’ are shown. It is well known that 
diffraction means a spreading out of the diffracted wavefront. The 
waveform of Fig. 4-3 diffracted and propagated out some distance might 
look like Fig. 4-4. Now at an even larger propagation distance the 
spreading out might be so severe that the wavefronts due to the object 
of interest and those due to the postulated ghost objects to either side 
of it will begin overlapping and interfering constructively or destruc- 
tively as in Fig. 4-5. Naturally the waveform of Fig. 4-5 looks quite 
different towards the boundaries of the sampling period from the one 
that would be obtained without any ghost objects and in fact is obtained 
and observed in a physically existing system as in Fig. 4-6. The observed 
waveform of Fig. 4-6 sampled and fed to the computer to be propagated 
back the correct distance using a DFT routine cannot reproduce exactly 
the central of Fig. 4-3, because of the limitations of the DFT routine. 

It is this fact that gives rise to the limitation in propagation distance 
of the transfer function approach to diffraction that was mentioned in 
the last section. There is an easy wasy to overcome this limit by 
increasing the guard space between the desired object and its ghost 
images by increasing the array size, if the implied increase in core 
memory needed and associated computing time can be afforded. In this 
investigation it was decided that a 64 by 64 element matrix was as 
large as could be afforded without running into time and memory problems - 
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Figure 4-4 

Waveform of Fig. 4-3 Propagated a Small Distance 
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Figure 4-6 

Central Object of Fig. 4-3 propagated the same distance as in Fig. 4-5 (no distortion) 



On this basis the largest useful propagation distance was found to be 
about 200 wavelengths for an object aperture with a linear dimension of 
the order of 15X. 

In image processing it is often desirable to have both the center 
of the input object and the zero frequency component (or undiffracted 
or d.c. component) of the spatial frequency spectrum centered in the 
array. When the DFT relationship of Eqs. (4-3) and (4-4) is used the 
zero frequency component appears in one comer of the transformed array. 
Andrews [Ref. 13] recommends multiplying the array to be transformed by 
(_l) n+ra before transformation, which will indeed shift the spectrum in 
such a way that its origin will be centered in the transform domain 
array. In this investigation a different approach was used because 
Andrews’ method still leaves an unwanted phase shift for an object that 
is centered in the middle of the array: the DFT routine asks for an 
object that is centered at the origin for an output without phase shift. 
The center of the array was defined to be at the point (l+N/2, l+N/2). 

A subroutine (called SHUFL, see Appendix) was written which interchanges 
the four quadrants of the array such that the array point (l+N/2, l+N/2) 
after rearranging occupies the point (1,1). (See Fig. 4-7) 

In interchanging the four quadrants of the array by a simple 
translation of the object, advantage was taken of the otherwise undesira- 
ble postulated existence of the ghost images. After transforming the 
array, the transform was once again operated on by SHUFL, which produced 
the desired effect of moving the origin of the angular spectrum to the 
center of the array (as previously defined) without any unwanted 
associated phase shift. 
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Array Before and After Use of Subroutine SHUFL 



B. FAST FOURIER TRANSFORM 



The discrete Fourier transform requires N 2 complex multiplications 
and additions for a one-dimensional record of N samples; this can be 
prohibitive for even moderately large N. The Fast Fourier Transform 
(FFT) is just a different way of computing a DFT (see {Ref. 14], but as 
it reduces the number of complex multiplications and additions on the 
order of N*log 2 N its importance can hardly be overestimated. The 
reduction in number of computations is even more dramatic for arrays 
of higher dimension. Reference [15] lists a few examples. Without the 
economy of the Fast Fourier Transform idea the DFT would be a useful 
tool, but to be used very sparingly because of the large time penalty 
to be paid. Since Cooley and Tukey published their report [Ref. 16] 
digital signal processing has received new impetus because transformation 
operations are now possible at a tolerable even though still not a very 
economical time expenditure. The FFT algorithm used in this investiga- 
tion (subroutine F0UR2) [Ref. 17], is included in the Appendix. The 
time needed for a transformation of a 64 by 64 element matrix using 
F0UR2 was less than 2 seconds. 
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V. THE DETECTION PROBLEM 



To recapitulate, the problem to be solved was to reconstruct the 
original object from its observed sound diffraction pattern without 
previous knowledge of the distance of the object from the observation plane. 

A. GENERAL OBJECT AND UNIVERAL RECONSTRUCTION CRITERION 

For an unknown object distance a given diffraction pattern can be 
due to an infinite number of ’objects'! It is therefore not possible 
without a prior knowledge to reconstruct a general object (not 
constrained in any way in amplitude and phase) by using a simple 
criterion function, i.e. one which does not contain the object 
distance in some coded form. The proof of this statement follows. 

Consider an object belonging to a rather general class U 
of objects, all of which are related to each other by some common 
amplitude and phase characteristics. Suppose also there is a 
reconstruction criterion W which when applied to all diffraction 
patterns A from all objects belonging to U successfully reconstructs 
the respective object U^. It is hoped that W is the 'universal 
reconstruction criterion'. 

Suppose a diffraction pattern A-^ appears at a distance z from U-^. 

Now take the diffraction pattern A* which appears at a distance z* in 
front of to be a new object to be reconstructed. A* gives rise 
to A^ at a distance z^-z* in front of A*. Application of W on A^ 
will always reconstruct U^ (according to the second hypothesis). 

Therefore it will never be possible to extract A* from A-^ by using 
the criterion W. So U cannot be the desired universal reconstruction 
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criterion. Since U and W were chosen arbitrarily in this example, 
one can generalize that no single simple criterion can be used to 
recover all objects. 

The expression ’simple criterion’ has been used in the above 
discussion for the reason that it was proposed to make two diffraction 
patterns of the same object using different wavelengths. Then, on 
backward propagation both images should coincide at the object location 
but not at different propagation distances, at least over a considerable 
portion of the array. This idea was not investigated any further 
because it was thought to circumvent the issue of the problem. The 
difficulties in this approach are obvious at once: The need of 

recording two diffraction patterns with the associated long time 
expenditure, the need of almost twice as much computer storage and 
time etc. Moreover it was hoped that a single diffraction pattern would 
contain enough information to find the original object if the class of 
objects to be reconstructed were somewhat limited. 

Recently a reconstruction time has been proposed by Kalra and 
Rodgers [Ref. 18] that works reliably on a class of objects, among them 
the point scatterer or source, which have a magnitude of the acoustical 
disturbance which reaches its maximum at the object location. Unfortunately 
this very simple criterion will not reconstruct successfully objects 
whose phase distribution departs much from a linear one. In the above 
reference [Ref. 18] the possibility of extending the analysis from a point 
source to a general object is suggested. Of course any object can be 
thought of as a (possibly infinite) sum of point sources and the 
diffraction pattern due to that object as resulting from the appropriate 
accumulation of point sources if the phase relationship between all 
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point sources is accounted for. But the assumption that maximum 
amplitude is a general reconstruction criterion can be disproven 
(without going to such exotic ’objects 1 as other diffraction patterns) 
by considering an uniformly insonified object, whose phase varies in 
a quadratic fashion such that it is advanced towards the edges and 
retarded in the center, i.e., a lens-type phase distribution. Then the 
sonic energy in propagation will concentrate near the focal point of 
this ’lens'. As the energy density there is much higher than at the 
uniformly insonified object (see [Ref. 19]), the amplitude of the sound 
field must be higher than at the object. This is due to the constructive 
interference introduced by the chosen phase relationship. The situation 
as above was simulated on the computer: A 16A by 8A slit with a phase 

distribution that had an equivalent focal length of 10A was insonified 
by a unity amplitude plane wave. The computed amplitude at the focal 
point was found to be more than a factor of 8 higher than the amplitude 
at the slit plane. So the proposed 'unbiased estimator' for the true 
object distance z* of above reference [Ref. 18] will always choose the focal 
plane as object location. It may be argued that a lens-type phase 
relationship is very uncommon in nature, but random phase distributions 
(diffuse scattering) are widely observed, and random phase relationships 
between the point sources that make up an object can also give rise to 
amplitudes of the disturbance which are higher at other planes, than 
at the object location. As example a 16\ by 8A unity amplitude object 
was taken to have random phase distribution between array points, 
uniformly distributed from 0 to 2 tt . The highest amplitude was found 
to be 2.1 times the value at the slit at a location 5A in front of the 
object. This must not necessarily be so; for other distributions the 
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amplitude may be highest at the object location (the linear distribution 
might be considered a very special case of a random distribution, as 
might be the parabolic distribution) , but it introduces an element of 
uncertainty into the reconstruction scheme which may well be unbearable. 

In the initial stages of this investigation it was hoped that cross- 
correlation of the diffraction pattern with a point source at the 
correct distance would give the necessary high output to discriminate 
between the object giving rise to the diffraction pattern and some other 
diffraction pattern due to this object. 

Suppose there is a unity amplitude point. source at the orgin. Then 
Up (x , y , o ) = 6(x,y,o) and 
Ap (fx 5 ^y 5 °) ” ^ 

If the disturbance due to this point source is propagated a distance z 
then 

Ap Z ( f x» f y» z ) = A p ( f x> f y > °) * ex P ^1- (Xf x ) (Xf y ) z ] 

= 1 • exp [j^-z /l-CAf^-CXfy)^] 

Now cross-correlation can be accomplished by multiplying the spectrum of 
one function by the complex conjugate of the spectrum of the other 
function and inverse Fourier transforming the result. 

Let the Fourier transform of the diffraction pattern ^(x^y >z) 
be given by A^(f^, z ) . Then the result of the proposed crosscorre- 

lation can be written as: 

R (x,y,o) * F -1 { A z (f x ,f y>z ) • exp[-j-^-z A- (Xf x ) 2 - (Xf y ) 7 ]} 
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which is recognized as the result obtained by backward propagation of 
the diffracted wavefront by a distance z, namely the original object. 

This method of maximizing the crosscorrelation cannot be used as it 
has been shown that the amplitude of the wavefront at the object need 
not be higher than everywhere else in its propagation field. 

B. OBJECTS WITH BOUNDARIES AND THEIR RECONSTRUCTION 

In the search for a class of objects which might be more useful 
than the one reconstructed by Kalra and Rodgers [Ref. 18] a brief study of 
psychopictorics , that subfield of psychophysics which deals with 
natural images as pictorial stimuli (see Lipkin and Rosenfeld [Ref. 20]) 
showed that the three most important psychopictorial variables are: 
contrast and border; shape and geometry; texture. Of these contrast 
and border is easiest to relate to what is actually available in a 
digitized picture. To treat shape and geometry would mean going to a 
pattern recognition system, probably by using a crosscorrelation scheme, 
with the attendant necessity of storing a large number of objects to 
be recognized. Trying to base a detection system on the texture of 
an object compared to that of its surroundings would probably be 
bound to fail with today’s still rather crude acoustical holography 
systems. Thus, trying to detect objects which are characterized by 
contrast and boundary seems most promising. In fact, many objects 
that are of interest to acoustical holography fall under this heading 
In that they are characterized by emitting or reflecting large amounts 
of sonic energy and being separated from the non-emitting (or less 
emitting) background by a rather sharp boundary or vice versa. 

The ’picture’ of a mine lieing on the ocean floor will show a sharp 
step in Intensity as will that of a bone imbedded in soft tissue 
and there are more examples that might be cited. 
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1. Objects of Interest 



It is hoped then that this class of objects will be general 
enough to warrant investigation. Therefore in all the following 
discussion it will be assumed that the object plane to be reconstructed 
from its diffraction pattern has associated with it (or ideally is 
characterized by) a sharp change in the emitted amplitude (ideally a 
discontinuity) separating regions of large amplitude from those of 
low amplitude. This change in amplitude will not be assumed to be 
confined to a single point or to a number of disconnected points 
because the word ’boundary 1 or ’edge 1 implies connection between quite 
a large number of points. In this sense then a point source is not 
a valid object for this discussion. A line object cannot be hoped 
to be detected, because it really has not any area associated with it; 
the same objection will hold for an object that is checkerboard-shaped 
with the amplitude changing from high to low and back again from one sampling 
point to the next one. The smallest object that one could hope to 
detect is composed of four sampling points of high intensity next 
to each other, arranged in a square. This probably would not be 
a great obstacle to th6 use of the proposed routine as it implies detection 
of at least all one wavelength by one wavelength objects if the Nyquist 
sampling rate is used. 

2 . Reconstruction Schemes 

The first reconstruction scheme for the given class of objects 
was very simple: As the object was assumed to have a sharp step in 
amplitude it was thought that the difference in amplitude between one 
sampling point and its neighboring sampling point would have a maximum 
at the object location and not in any other plane. Therefore that 
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plane which contains this maximum would be taken to be the object 
plane and the reconstructed diffraction pattern in this plane would be the 
image of the original obj ect • Unfortunately the above assumption was 
proven to be not valid for the Nyquist sampling rate if the phase 
distribution over the object departs markedly from a linear one, by 
assuming such an object and simulating its diffraction on the computer. 
Random phase objects as well as such with a lens-type phase relationship 
showed maximum differences between some sampling points which were 
higher (even very much higher in the case of the parabolic phase object) 
in other planes than the discontinuity in the object plane. As examples: 

A 16X by 8X unity-amplitude object with uniformly distributed random 
phase had a maximum difference of 1.73 at a distance 3 wavelengths in front 
of the object. The same object with lens-type phase distribution and 
a focal length of 10X showed a maximum difference of 3.85 in the focal 
plane. At first this was not taken to invalidate the whole idea, 
because theoretically by finer and finer sampling it would be possible 
to reduce the maximum difference between sampling points in any other 
plane so far that it would be smaller than the given discontinuity in 
the object plane, which would not be influenced by finer sampling. 

But apart from the obvious difficulty in an actual system to sample 
fractions of a wavelength apart (the detector would have to have very 
small area indeed) the hoped for gain could only be realized at an immense 
increase in amount of data and additional computation time and core 
memory necessary to handle that many data. If one samples finer without 
increasing the number of samples (which might be possible for very 
small objects) nothing is gained. This hinges on the fact that by 
sampling above the Nyquiest rate too few frequency components will be 
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available to faithfully reconstruct the original object. (see Sect. IV) 

The image will lose its clear amplitude step and instead will develop 
a lower amplitude slope and overshoots; in this way the characteristic 
that discriminates between actual object and its diffraction pattern is 
lost. Some other way has then to be found to detect the edge of the 
intensity distribution. 

In the above scheme no use was made of the fact that an edge 
consists of more than one spot where the amplitude changes rapidly, 
that a number of such changes have to occur at adjacent sampling 
points and that these places of rapid change have to separate a 
region of high amplitude from one of low amplitude. 

The problem of enhancing or detecting an edge in a digital picture 
has been treated in the literature, e.g., by Rosenfeld, Lee, and 
Thomas in [Ref. 21]. The usual treatment is to take differences of averages, 
e.g. for a one-dimensional case 



dmk = — 
m 



ak+m + •••+ ak+1 - ak -••• -ak_ m+ i 

where a^ is the amplitude of the point being investigated. The 
difference d^=a^ + ^~a^ is sensitive to every local fluctuation, 
whether it occurs at an edge or not. The difference d^ for large m 
has the disadvantage of detecting even sharp transitions in a blurred 
manner since it registers a difference for all k such that some of its 
terms are in one segment, some in the other segment. A simple method 
of overcoming these difficulties is using a product of d^’s, which 
Involves both small and large m’s. This product can only be large 
if all its factors are large. For the two-dimensional case this edge 
detection scheme naturally becomes more involved. Operations of the form 
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^rs 5 nk r (2s+l! 



\ ^ a h+r 5 k+s 



+ • • *+a 



h+r ’k-s 



+• • *+a. 



h+1 k+s 



+• • • 




h-r+l 5 k+s 



+ ^-r+l’k-s^ I 



are performed and products of the results for various large and small 
r’s and s's are formed. Again the product can only be relatively 
large if all its factors are large. The example d shown above 

is used for the detection of a horizontal edge, the parameter s increas- 
ing results in the sensitivity towards orientations departing from the 
horizontal being accentuated. Each possible edge orientation has 
to be treated by a different combination of differences of averages, 
obviously a time-consuming task for an array of say 64 by 64 sample 
points. The edge detection operation would normally then be followed 
by a curve detection routine to filter out any spurious, output, as the 
output of the edge detection step will be noisy for an edge which is not 
perfectly smooth and straight. The curve detection idea is based on the 
fact that the values resulting from the edge detection should be high 
on moving along the edge and low in a direction perpendicular to it. 

In other words a high point on an edge should have neighboring high points 
In approximately opposite directions. Points which do not satisfy this 
criterion are regarded as spurious and discarded. For a more profound- 
discussion of this topic reference is made to Rosenfeld, Thomas and 
Lee [Ref. 22]. Again, the time needed for this step in the recognition 
process cannot be regarded as negligible, because it involves checking the 
neighborhood of N 2 points and a quite involved decision-making process. 

To receive some information on the computation time needed using 
the above method of edge detection a first approximation to it was 
developed and tested. In order to make the operation d rg , n ^ relatively 



62 



insensitive to edge orientation s was taken to be 1. Differences 

d 21’hk ,d 4l’hk and d 81’hk were forTned and multiplied together 
for the horizontal and vertical edge orientation and the products 
for the horizontal and vertical operation added together. It was 
hoped in this way to overcome the necessity of investigating all 
different orientations. The result obtained with this technique was 
promising. (The test was performed using a circle, whose area was 
filled with high values, as a boundary comprising all edge orientations 
as input to the routine. Input and output are shown in Fig. 5-1 and 
Fig. 5-2). Nevertheless this method had to be shelved because of the long 
computation time needed, over 20 seconds for each 64 by 64 element 
array. To this time would have to be added the time needed for the 
previously discussed curve detection algorithm also. 

The difficulty of having to treat each orientation separately 
is analogous to that which is encountered in trying to accomplish edge 
detection by a cross-correlation scheme, which also is very sensitive 
to a relative rotation between the reference object and the object to 
be tested, and additionally has to account for different phase relation- 
ships also. 

In order to surmount these difficulties due to various edge orienta- 
tions a new, rather simple (and thus more error-prone) procedure was 
worked out. In this algorithm first the largest difference between any 
two adjacent array points in the plane under investigation is found. 

This is used to set two thresholds, an upper (THRUP) and a lower one 
(THRLO) . These thresholds are kept variable, so that they can be adjusted 
according to the kind of object one is looking for (if such previous 
information is available). A setting of 60% of the maximum difference 
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Figure 5-1 

Amplitude Contour Plot of Input to 
Edge Detection Scheme to be Tested 
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Figure 5-2 

Amplitude Contour Plot of Output of 
Edge Detection Scheme to be Tested 
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between array apoints for the upper threshold and 20% for the lower 
one seems in general to give good results. An edge in one direction 
is assumed to exist at an; array point (i,k) if for a downgoing step 

^ a i *k‘" a i»k+l^ is S reater th an the upper threshold and both I a -j_ >k+l” a i *k+2 I 

and |a., -a . , | are less than the lower threshold. This criterion is 

1 *k-l i *k 

analogous to the edge detection routine described earlier with S =* 0 and 
the product operation replaced by its logical equivalent f AND f : The 

output in both cases will be positive ( T YES T in this case, large in the 
former) if all the inputs (differences) to the operation are positive 
('YES* or large). 

This edge criterion is used twice at each array point, for horizontal 
and for vertical differences. In this manner corners as points belonging 
to two edges at the same time are weighted more than simple edge points. 
(The deemphasis of comers is one unfortunate feature of the earlier edge 
detection scheme). Also, a distinction is made between upward steps 
^ a i’k+l ^ a i J k^ must g reater than the upper threshold, while the other 
two conditions stay unchanged in this case) and downward steps. This 
feature is used in that part of this scheme which assigns different values 
to each edge point according to whether it is connected to other edge 
points or not: The step heights for an array point (i,k) are computed in 

both the horizontal and vertical directions; values of positive going 
steps are kept separate from negative steps; then the four array points 
next to (i,k) which have been treated previously in the sequential edge 
search are tested for the highest absolute value of assigned 1 edgemeasure f . 
If this is due to an upward going edge the height of the positive step at 
the point in question is added to it and stored at (i,k). The treatment 
is analogous if the highest value stored next to (i,k) has been found 
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from a negative edge; the negative step height at (i,k) is added to it 
and kept at (i,k). If the highest adjoining edge measure is positive, 
but only negative steps occur at (i,k) this negative step height is 
stored at (i,k) without any addition of previously computed edge measure. 
The same thing happens if the adjoining array points contain zero. 

Zero is stored at (i,k) regardless of adjoining edge measures if the 
situation at (i,k) and its vicinity does not meet the edge criterion 
using the upper and lower thresholds. For a more exhaustive treatment 
of all different cases that may arise see the computer flow diagram for 
subroutine ADERIV in the Appendix. After the whole array has been 
treated in this manner, all edge measures computed are added in absolute 
value; the result is taken to be the overall edge measure of the given 
array. 

It must be admitted that the only justification for this scheme 

that can be given, is that it seems natural in looking for edges to 

discriminate between large and small differences between adjoining 

points and that longer connected boundaries should be weighted more than 

short ones: three points of rapid change next to each other might be 

considered at least as a beginning of an edge, whereas two such points 

could be just a random occurrence. The normal height of the thresholds 

as given is found entirely from practice with this routine and has no 
* 

profound mathematical background. It may be argued that in the 
decribed algorithm too much importance is attributed to connection 
between edge points and too little to the height of the single, steps . 

But it has been observed that if an array containing fewer, large steps 
is compared with another array with many smaller steps next to each 
other, the latter is more likely to b e the looked for object than the 
former. Besides that, a shift in the relative weight balance between 
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’connection’ and ’step height’ is easily accomplished by multiplying 
the largest absolute value of edge measure adjoining the point (i,k) by 
some factor less than one before the addition and storage at(i,k). 

The algorithm can easily be made more discriminating against steps 
which do not belong to an actual edge by taking into account not merely 
two differences next to the point in question but investigating more 
points in its neighborhood. Then it must be borne in mind that as more 
small differences in one direction are required to satisfy the step 
criterion the smallest detectable object will become larger. Another 
means of increasing the rejection capability is by increasing the 
upper threshold and reducing the lower threshold if that is compatible 
with the kind of object to be investigated: setting a lower threshold 
of 0.1 of the maximum difference, if the amplitude of the background of 
the object is uniformly random distributed between zero and 20% of the 
object discontinuity, could reject quite a number of bona fide edge points 
thus reducing the probability of detection. It is inadvisable to work 
with an edge detection routine which is too discriminating from another 
point of view also: this routine might not pass the diffraction patterns 
a small distance behind and in front of the object plane with a high 
enough edge measure to make them obvious. Thus in stepping through the 
volume of interest the computer would have to investigate very near to 
the actual object to receive an indication of its presence or it would 
be missed altogether and the computer would settle for some other than 
the true object distance and output the diffraction pattern there as 
the reconstructed object. This would mean that many small steps would 
have to be taken, which implies very long computation times. This case 
could not appear if the overall edge measure as computed above were 
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steadily increasing towards the actual object location. Then it would 
be just a matter of finding the maximum edge measure and the object 
would be detected. Unfortunately, the situation is not that simple. 
Figure 5-3 shows the computed edge measures for a 16 by 8 wavelength, 
unity amplitude object with random phase, versus propagation distance 
in steps of 2.5 wavelengths. It shows quite a number of planes in 
which the edge measure has a relative maximum. This is not only due 
to the way in which the edge measure is computed but also due to the 
rapid variation of the diffracted wavefront with distance. 

A reconstruction method then has to s tart from a diffraction pattern 
and taking advantage of all previous information (which might have been 
obtained by rangegating) propagate this diffracted wavefront back to some 
starting distance. From there to some given final distance the wavefront 
will be propagated back in as small a number o f steps (as large a step 
size) as possible without missing the object altogether, meaning that 
the step size has to be so small that an indication of a possible object 
is received by a relative edgemeasure maximum less than half a step size 
behind or in front of the actual object location. In the proposed 
algorithm this is accomplished by setting the parameters ZFIRST for the 
starting distance, ZLAST for the final distance and NRSTEP for the 
number of steps to be used. The choice of step size for this first pass 
through the volume of interest is thus left to the operator. For an 
object of 16 by 8 wavelengths and Nyquist rate sampling the step size 
should not be chosen larger than about 2 wavelengths. Larger objects 
can bear a larger step size. 



69 




Figure 5-3 

Edge Measure of Diffraction Pattern of Random 
Phase Slit Versus Propagation Distance 
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All relative maxima in overall edge measure observed on this first 
pass have to be stored for future reference, together with the propaga- 
tion distances at which they occur. 

Different treatments of the result of this first pass are possible 
now. For instance the reconstructed wavefront at all locations of a 
relative maximum might be displayed on a plotter, or graphics terminal 
and the decision as to the actual location(s) of object(s) left to the 
experience of the operator, who would then investigate those locations 
closer. The method chosen for the proposed algorithm differs from the 
above by making further use of the computer's decision-making 
capabilities. 

After the first pass through the whole distance of interest the 
vector of edge measure maxima is reordered according to their size and 
those possible object locations discarded, whose edge measure is only 10% 

t 

or less of the maximum edge measure encountered on the first pass, it 
being assumed that closer investigation of these locations will not 
yield the looked for object. The remaining locations are investigated 
in a fixed step size which is set equal to the sampling increment, i.e, 
to one half a wavelength if the Nyquist rate is used, from a distance 
one half of the original step size in front of the relative maximum to 
one half that step size behind this location. The absolute maximum in 
edge measure in this interval is found and stored in a different array. 
After all relative maxima resulting from the first pass have been 
investigated the array of second pass maxima is scanned for the overall 
maximum and those locations discarded where the edge value is below 
a certain fraction of the overall maximum. This is accomplished by 
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setting a threshold THRCON to a fraction of 1.0: one would not expect 
maximum of less than 0.2 to 0.3 of the overall maximum to be due to 
actual objects. Also incorporated into the algorithm was the possibi- 
lity of investigating fewer relative maxima resulting from the first 
pass by reading in a value for NRMAX, e.g., if it is known that the 
object exists only in a single plane parallel to the observation plane. 
This option should be used with caution because there is no guarantee 
that the highest relative maximum after the first pass will be associated 
with the actual object for all objects (although this is true for many 
simple objects). 

The results of propagating the wavefront to the locations of the 
remaining second pass maxima are then plotted on the CALCOMP - plotter, 
either in form of a contour plot of the computed amplitude or as a 
surface composed of the values of amplitude at all array points, from 
which the hidden lines are removed (other graphical output devices can 
also be used). The contour plot will give a much better pictorial 
representation of the outline of the object, but a good part of the 
amplitude information is lost by plotting only a number of contour lines. 
The three-dimensional representation is t o be preferred on the basis 
of information contained in the plot (it also suffers from loss of 
information because it is not possible to see through the surface) 
but it takes much longer computation times and requires much more 
core memory than the contour plot. Also it must be remembered that 
the apparent dimensionality is not due to any three-dimensional object 
of that shape but rather represents the relative amplitude of the 
wavefront in a single plane. Subroutines used for this output are for 
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the contour plot CONTUR [Ref. 23] and PR01 for the surface plot 
[Ref. 24]. For a 64 by 64 element array to be treated CONTUR takes 
on the order of 50 sec if six contour lines are plotted; PR01 needs 
for the same array from 50 sec up to 10 min depending on how complicated 
the surface to be plotted is. 
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VI. TESTS AND RESULTS 



Part of the tests performed were designed to check the propagation 
algorithm, the greater part were used to try the automatic reconstruction 
scheme. For this latter purpose cumputer generated diffraction patterns 
and one set of actually observed data were used. 

A. TEST OF THE PROPAGATION ALGORITHM 

As the spatial frequency approach to diffraction is limited in the 
maximum useful propagation distance due to the presence of the postulated 
ghost images it is not possible to apply the usual test for diffraction 
algorithms in which the amplitude distribution at a very large progaga- 
tion distance has the functional form of a Fourier transform of the 
diffracting aperture. Instead the fact was used that illuminating an 
aperture with a spherical wave converging to a point P Q will give rise 
to an amplitude distribution in a plane perpendicular to the direction 
of propagation and containing P Q , which again is simply related to 
the Fourier transform of the diffracting aperture. [Ref. 7]. 

A 10X by 10X square aperture was simulated to be insonified by 
a spherical wave converging to a point 40X behind the aperture. Using 
the proposed propagation algorithm and propagating the diffracted 
wavefront 40X resulted in the output of Fig. 6-1. It is readily 
apparent that this output has the functional form of (sin(x) /x) • (sin(y) /y) 
(which is the Fourier transform of a uniformly insonified rectangular 
aperture). The fact that the surface of Fig. 6-1 is not very smooth is 
due to the relatively small number of surface points. Due to the 
presence of the ghost images the output is deteriorated towards the 
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boundary of the presented surface (the amplitude does not fall off 
as fast as it should in following a (sin(x) /x) • (sin(y) /y) functional 
relationship) . 

It is of interest to note that even after about 30 sequential 
propagation steps the result differs from that obtained in one step 
only in the fourth decimal place. This good agreement could never be 
expected if the Fresnel equation were used. 

B. TESTS OF THE RECONSTRUCTION SCHEME 

1. Tests Using Computer Generated Diffraction Patterns 

Most of the tests of the automatic reconstruction scheme were 
performed using computer generated diffraction patterns. They have 
the advantage of being perfect in the sense of having no distortion 
introduced by noise in the detection system, by reverberation and by 
scan non-uniformities. In all tests involving computer generated 
diffraction patterns the distance from object to observation plane 
was chosen arbitrarily (subject to the constraint of being small 
enough for the results of the spatial frequency approach to diffraction 
to be valid) . The numbers of steps given in the following discussions 
are those for the first pass through the volume of interest (i.e, 
those given by setting the input NRSTEP) . The number of steps for the 
closer investigation of relative maxima in edge measure is computed by 
the detection algorithm (as discussed in the last section) . All com- 
putation times are without those necessary for the output subroutines 
(CONTUR or PR01) . 

Case a: The simplest object whose reconstruction was attempted was a 

16A by 8A uniformly insonified rectangular slit with all its 

sample points in phase thus simulating a specular object. The diffraction 
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pattern was observed 175. 3X in front of the object; in the recon- 
struction the closest distance to be investigated was set to 160X, 
the furthest distance was 192X; this volume of interest was investi- 
gated in 12 steps. A reconstruction (see Fig. 6-2) was obtained at 
a distance of 175.25 wavelengths from the observation plane, 
i.e. only 0.05X away form the original object. The time needed was 
about 2.6 min. 

Case b: The same object as in case a was assumed with the additional 

difficulty of a background insonification with an amplitude which was 
randomly distributed between zero and 0.2 of the object amplitude. 

Also the phase of the background insonf ication was taken to be 
uniformly random distributed between zero and 2tt, The distance 
between object and observed diffraction pattern was 177. 2X. The 
volume of interest and number of steps were as in case a. The obtained 
reconstruction of Fig. 6-3 was located 0.05X behind the actual 
object location. In order to show the random amplitude variation of the 
background here the output is given in form of a 3-dimensional surface. 
Computation time was about 2.6 min. 

Case c: An object of the same size and shape as in case b, was treated 

in this case. The random background, too, was the same as in case b. 
This case differed from case b, in the assumption of a parabolic phase 
relationships for the object which gave it an equivalent focal length 
of 10 wavelengths. The observation plane was located 177. IX in front 
of the object. Initial and final distances remained unchanged from cases 
a and b, as did the number of steps. An object reconstruction was 
obtained at a distance of 177.25, requiring a computation time of 3.3 
min. 
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Figure 6-2 

Amplitude Contour Plot of 
Reconstructed Specular Object 
(Case a) 
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Case d: To show that the proposed reconstruction scheme works also for 

other than rectangular objects a shape was chosen for the input object 
which can best be visualized by viewing the reconstructed output of 
Fig. 6-4. Background and phase distribution of the object are the same 
as in case c. Propagation distance to the observation plane was 177. IX. 
Initial and final distances for the investigation were set to 160.0 and 
190. .0 wavelengths respectively, the investigation was performed in 12 
steps. The reconstructed object (see Fig. 6-4) was located 177. 25X 
from the observation plane. Time needed for the reconstruction was 
about 4 min. Figure 6-5 shows the same output as Fig. 6-4, only in the 
3-dimensional surface plot. This figure was included here to illustrate 
some of the shortcomings of this output scheme: Even for a relatively 

simple shape as that from Fig. 6-4 the 3-dimensional plot does not 
show the object shape clearly. Moreover it exhibits certain errors 
due to the rapid variation of the amplitude function which cannot 
sufficiently be accounted for in the interpolation scheme. The contour 
plot (Fig. 6-4) on the other hand, does not show the random variation 
of amplitude in the background insonif ication while exhibiting quite 
clearly the shape of the object. 

Case e: For this case a 16X by 8A rectangular uniformly emitting 

object was assumed to have a uniformly distributed random phase distri- 
bution between zero and 2 tt whcih simulates a diffuse object. No back- 
ground sound was taken to exist. The diffraction pattern was calculated 
for a observation distance of 177. 5X. The volume of interest remained 
unchanged between 160X and 190X, as did the number of steps in the 
investigation. One would expect this case to be harder to solve due 
to the strong deterioration of the best possible reconstructed object 
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Figure 6-4 

Amplitude Contour Plot of Reconstructed Parabolic Phase 
Object with Random Phase, Random Amplitude Background (Case d) 
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Figure 6-5 

Surface Plot of Amplitude of Reconstructed Parabolic Phase 
Object with Random Phase, Random Amplitude Background (Case d) 



(see Fig. 3-6). The computer found two possible object locations. 

Plots for the diffraction patterns at these locations (177. 5 A and 
188. 25A) are shown in Fig. 6-6 and 6-7 respectively. The first of these 
locations is the actual object location and the plot Fig. 6-6 the 
best output that can be expected. Computation time was 5.5 min. It 
is surprising that the edge-detection routine does not discard the 
array of Fig. 6-7, but instead assigns it a value of overall edge- 
measure which is 70% of that of the actual object. Here it would be 
necessary for a human operator to use his pattern-recognition capability 
to discard the second object hypothesis. 

Case f : As the last test for the reconstruction of objects which exist 

only in a single plane a 16A by 8A object was simulated with a 
background insonif ication whose amplitude was uniformly random distributed 
between zero and 0.2 and whose phase was assumed to vary uniformly between 
zero and 27 t. The object itself had an amplitude which consisted of 
a constant term of 0.9 to which was added a randomly varying value 
which was uniformly distributed between zero and 0.2. Also, a random 
phase variation between zero and 2 tt was assumed for the object. The actual 
distance between object and observation plane was 177. 5X. The initial 
distance, final distance and number of steps were left the same as in 
case c. The computer required 5 min. to calculate two tentative object 
locations at 177. 5A and 163. 75A. The diffraction patterns at these 
locations are shown in Figures 6-8 and 6-9 respectively. The edge- 
measure assigned by the computer to the spurious output of Fig. 6-9 was 
0.75 of that of the correct output. 

From all these tests it would seem that one could expect the 
computer to be able to reconstruct any object which has a sharp edge 
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Figure 6-6 

Amplitude Contour Plot of Reconstructed Random Phase Object 

(Case e) 
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Figure 6-7 

Amplitude Contour Plot of Spurious Output 
(Case e) 
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Figure 6-8 

Amplitude Contour Plot of Reconstructed Random Phase Object 
with Random Amplitude Variation. Background Assumed to have 
Random Phase, Random Amplitude (Case f) 
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Figure 6-9 

Amplitude Contour Plot of Spurious Output 

(Case f) 
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and exists only in a single plane with rather high probability. 

Naturally each output would have to be scrutinized by an operator 
such that spurious output which could appear in addition to the actual 
object reconstruction would be discarded. The proposed edge detection 
routine is not infallible, but it performs 1 its operation in rather 
short time which was the main reason for setting it up in the existing 
manner. 

The next set of tests performed consisted of trying to reconstruct 
input objects existing in two different parallel planes. As in backward 
propagation the real image is reconstructed, the first difficulty is 
immediately obvious: the diffraction pattern converges to that existent 

in the plane of the object on backward propagation, but the sound waves 
emanating from the object do not stop on the reconstructed object but 
are propagated farther back on the next propagation step. This poses 
no great problem for an object that exists only in a single plane parallel 
to the observation plane. But consider an object consisting of two 
extended sound sources in two different parallel planes. The sound from 
the source farther from the observation plane is propagated into the 
plane of the nearer source there to be absorbed, phase shifted and/or 
diffracted, as the case may be. The sound field that then exists in the 
plane of the front source is propagated to the observation plane. On 
backward propagation the field converges to that previously existent 
in the front plane. So far, no fault can be found with the described 
situation. But on further backward propagation (the second source should 
also be reconstructed) the sound field introduced by the front source 
is also propagated back into the plane of the farther object, there to 
give rise to a disturbance which did not exist in actuality. In fact 
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this is not too different from the case of the optical hologram and the 
eye because the eye converts the virtual image to a real image on the 
retina. By focusing on a part of the image that is in a plane farther 
away, that part which is nearer is still seen, in a blurred fashion, 
in that same plane, but the brain rules out the hypothesis of a 
composite sharp and blurred object: The unwanted information is largely 

ignored. How recognition is accomplished, i.e., why certain object 
hypotheses are accepted, other not, remains unknown to date, even 
though the necessity of experience is generally assumed (apart from 
the obvious advantage of stereoscopic vision if the distance is not too 
large). Huang [Ref. 25] and Ichioka, Izumi and Suzuki [Ref. 26] show a 
method of generating diffraction patterns due to such objects. The 
disturbance due to the back object is propagated into the plane of the 
front object. There that part of the resulting disturbance that 
corresponds to the extension of the front object is replaced by the 
disturbance due to the front object and the resulting composite 
disturbance is propagated into the observation region. This scheme 
was used for generating the necessary diffraction patterns in this 
investigation. 

Case g: The simplest object of this category consisted of two 16A 

by 8A rectangles, assumed to be specularly emitting (i.e,, all sample 
points were in phase). These objects were located in planes which were 
15 \ behind one another. One rectangle was rotated 90° with respect to 
the other around the propagation axis and displaced in such a manner 
that the projections of both rectangles along the axis of propagation 
overlapped by only a small amount (about 0.1 of the area of one 
rectangle). The actual object configuration can very well be 
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visualized by studying Figs. 6-10 and 6-11, which are the results of 
the computer reconstruction. The diffraction pattern was observed 
100.4 wavelengths in front of the front rectangle. The volume of 
interest was taken to extend from 85X to 130X from the diffraction 
pattern. Therefore the number of steps in the investigation had to be 
increased to 30. The reconstructions obtained were located at 100. 25X 
and 115. 25X from the observation plane, i.e., 0.15X in front of the 
first and second rectangle respectively. The total computation time 
was about 8.5 min. 

Case h: The same geometry as for case g was assumed but here both 

rectangles had a random phase distribution uniformly distributed from 
zero to 27 t. The volume of interest and number of steps for the 
investigation were not changed from case g, but the observation plane 
was located 99.63X in front of the front rectangle. The computer 
obtained 8 object locations (at 93.53X, 99.69X, 103. 85X, 106. 95X, 

116. 26X, 120. 91X, 123. 79X and 126. 12X), a clear indication that this 
reconstruction was not quite successful. The array corresponding to 
the actual location of the front rectangle, Fig. 6-12, was only third in 
overall edge-measure. An output which is typical for the spurious 
objects computed is shown in Fig. 6-13. The amplitude distribution of 
Fig. 6-13 was the fourth highest in edgemeasure; it is computed at a 
distance of 7.32X behind the front object. Figure 6-14 shows the output 
which is the lowest of all eight computed in edgemeasure. Its generating 
diffraction pattern is located 1.6X behind the second rectangle. It is 
not evident that this output is strongly related to the second rectangle 
(a comparison with Fig. 6-11 is helpful). The existence of the front 
rectangle can be seen quite clearly in Fig. 6-12 (compare with Fig. 6-10). 
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Figure 6-10 

Amplitude Contour Plot of Reconstructed Diffraction 
Pattern in Front Plane of Double Specular Object 

(Case g) 
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Figure 6-11 

Amplitude Contour Plot of Reconstructed Diffraction 
Pattern in Back Plane of Double Specular Object 
(Case g) 
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Figure 6-12 

Amplitude Contour Plot of Reconstructed Diffraction 
Pattern in Front Plane of Double Random Phase Object 

(Case h) 
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Figure 6-13 

Amplitude Contour Plot of Spurious Output 

(Case h) 
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Figure 6-14 

Amplitude Contour Plot of Reconstructed Diffraction Pattern 
Near Back Plane of Double Random Phase Object 

(Case h) 
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On studying Fig, 6-12 and 6-13 it is not surprising that the computer 
cannot clearly recognize the output of Fig. 6-12 as the actual object. 

The deterioration of the output due to the loss of evanescent waves is 
very apparent in a situation as complicated as the one treated in this 
case. 

Case i: In order to investigate the effects of varying phases between 

both rectangular objects and those of larger amounts of overlap, two 
rectangles of unity emission amplitude were simulated. These were 
assumed to be emitting in phase, but their separation in the direction 
of propagation was taken to be an odd number of half wavelengths, 
namely 14. 5X. Both rectangles were of size 16A by 8A, one was rotated 
90° versus the other around the propagation direction and both were 
centered in the input array. The overlap thus became one half the area 
of one rectangle. The diffraction pattern was observed at 100. 4A from 
the front object. Investigation took place in 30 steps between 85X and 
130X from the observation plane. The computer generated two possible 
object locations at 100. 75X (Fig. 6-15) and at 113. 5X (Fig. 6-16) from 
the diffraction pattern. This took 6.2 min. of computation time. The 
front rectangle being reconstructed 0.35X behind its actual location 
shows up very well. Due to the destructive interference between front 
and back object the location of the second rectangle could not be found 
exactly but the shape of the second object is quite apparent in Fig. 6-16. 
Removing the front object from the diffraction pattern, once it is 
recognized, can be helpful for finding the location of the second object 
in a case like this (see case j, below). 
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Figure 6-15 

Amplitude Contour Plot of Reconstructed Diffraction 
Pattern in Front Plane of Double Specular Object 

(Case i) 
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Figure 6-16 

Amplitude Contour Plot of Reconstructed Diffraction 
Pattern Near Back Plane of Double Specular Object 

(Case i) 
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Figure 6-17 

Amplitude Contour Plot of Reconstructed Diffraction 
Pattern in Back Plane of Double Specular Object 

(Case j) 
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Figure 6-18 

Amplitude Contour Plot of Diffraction Pattern in Front Plant of Double 
Specular Object, Obtained After Removing Back Object. (Case j) 
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Case j: A final test using computer generated diffraction patterns 

was performed with an input geometry as in case i, with the difference 
that the distance between rectangles was set to 15. OX and both were 
taken to be specularly emitting in phase. Distances to the observation 
plane and from these to the first and last plane of interest remained 
unchanged as did the number of investigation steps. The computer 
recovered the location only of the back rectangle at 115. 25X 
from the observation plane (see Fig. 6-17). This was accomplished 
using 7.5 rain of computer time. It is surprising that the plane of the 
nearer rectangle was not found. In a case like this it may be helpful 
to remove the clearly recognizable object from the output by setting 
a portion of the reconstructed array to zero and to use the remaining 
diffraction pattern to try to recover the rest of the composite 
object. This was done here: The clearly indicated back rectangle was 

zeroed out and the remaining disturbance propagated forward 25X. A new 
volume of interest was defined to extend from 5X in front of the 
new observation plane (as the diffraction pattern there was not 
displayed, it was not on actual observation plane) to 40X , which was 
investigated in 25 steps. This took 6.3 min. of computation. The 
computer recovered the front object 0.15A in front of its actual location. 
The resulting diffraction pattern there is shown in Fig. 6-18. The 
fact that the center portion of the rectangle is distorted is attributable 
to the removal of the center portion in the back object plane. 

The obtained results from this series of tests imply that automatic 
reconstruction of objects from their diffraction pattern is feasible. 

Of course for this to be accomplished not only the objects but also 
their reconstructions must meet the criterion which is applied to make 
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the object - no object decision. This is the reason why in the above 
tests certain objects could not be recognized by the computer. Neglecting 
the evanescent wave destroys a larger or smaller amount of edge 
information, depending on the object. This has to be kept in mind 
when judging the merits of the proposed scheme. 

2 . Testing the Reconstruction Scheme with Experimental Data 

As the input for the final test a set of data from an actual 
diffraction pattern was used. This diffraction pattern had been observed 
by C. Griggs (see [Ref. 27]) for a different purpose. The object giving 
rise to the diffraction pattern was a square transducer of 2 inch side 
length. The active area was 1.75" by 1.75". The pattern was sampled 
with a linear transducer approximately 73X away from the transducer 
by use of a scanning system. These data suffer from a few minor details 

Csee [Ref. 27])> the most important of these being that only one half of 

the diffraction pattern was observed to conserve time. It was therefore 
necessary to fill the rest of the input array by a mirror image of the 
actually observed data. This explains the unnatural symmetry that can 
be observed in inputand output (Fig 6-19 through 6-24). Another 
drawback of this set of data, which has to be taken into account, is the 
fact that the difference between adjacent sampling points is 0.9X, 
i.e., the sampling was performed at a rate very much below the Nyquist 
rate; aliasing effects are therefore to be expected, mostly at high 
spatial frequencies which means that distortion of edges by enhancement 
or deemphasis is probable. But the most important reservation against 
these data is not due to the manner in which they were recorded but 
rather to the object which generated this diffraction pattern: As 

the transducer had to be clamped at its boundary to get a watertight 
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seal it is quite possible that it does not constitute a good object for 
the proposed reconstruction scheme. One cannot assume that its 
amplitude dropped abruptly from a rather high value to a very low one, 
but would think that his transition would be more gradual. 

Nevertheless the test was performed with these data as they 
were readily available. 

Amplitude and phase of the diffraction pattern after reading 
in and reflecting about the center are shown in Figs. 6-19 and 6-20. 

The volume of interest was defined to be between 60A and 80A 
from the observed diffraction pattern. This valume was investigated in 
8 steps. This investigation took about 4.4 min. Two possible object 
locations were found at 74.66 \ (edge-measure 29.927) and at 67.68X 
(edge-measure 26.438) from the observation plane. At first sight there 
is nothing in the amplitude display (Fig. 6-21) from the location with 
the highest edge-measure, to prefer it as recovered object to that 
of Fig. 6-22, which shows the amplitude distribution at 67.68A from the 
observation plane (or for that matter from the input amplitude 
distribution), except for the fact that it appears close to the given 
nominal distance of 73X. In fact the amplitude pattern for such an 
object does not change very much in the volume of interest. The 
plausibility that the correct location was assigned the highest 
edge-measure becomes quite high when the phase distribution in Fig. 

6-23 (at74.66A) and in Fig. 6-24 (at 67.68A) are studied. It can be 
seen that in Fig. 6-23 all points on the face of the transducer have 
the same phase (or very nearly so), whereas in Fig. 6-24 a distinct 
variation can be observed. 
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Figure 6-19 

Amplitude Contour Plot of Experimental 
Data; Input to Reconstruction Routine 
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Figure 6-20 

Phase Contour Plot of Experimental Data; 
Input to Reconstruction Routine 
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Figure 6-21 

Amplitude Contour Plot of Reconstructed 
Diffraction Pattern with Highest Edgemeasure 
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Figure 6-22 

Amplitude Contour Plot of Second Obtained Reconstruction 



107 




Figure 6-23 

Phase Contour Plot of Reconstructed Diffraction* 
Pattern with Highest Edgemeasure 



108 




Figure 6-24 

Phase Contour Plot of Second Obtained Reconstruction 
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The result of this test seems rather inconclusive* but this is 
due more to the nature of data available rather than the proposed 
reconstruction scheme. 
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VII. CONCLUSION 



Computer reconstruction of objects from their diffraction patterns 
has certain advantages versus optical reconstruction, notably the 
avoidance of longitudinal distortion due to the difference in acoustical 
generating and optical reconstruction wavelength. The main drawback 
of computer reconstruction is the need of knowing exactly the distance 
of the object from its diffraction pattern. For an unknown object 
distance one would have to use many small steps of backward propagation, 
display the diffraction pattern at each location and have an operator 
make the decision as to the object location according to his preconceived 
ideas of what an object looks like. As the manipulation of data into 
form useable for output takes up considerable amounts of computer time 
it will be helpful if only a small number of possible objects is 
displayed. Automatic object recognition by computer can be an aid here. 

Automatic object detection by a computer routine can only succeed 
for objects belonging to that class for which the routine was written, 
since it is impossible to find a criterion which will reconstruct all 
objects. A particularly important class of objects is that which 
includes objects with sharp boundaries and edges, i.e. , rapid changes 
in amplitude separating regions of high emission from those with much 
lower amplitude. Therefore the reconstruction scheme presented in this 
report was tailored to detect edges and boundaries, without using 
great amounts of computer time and core memory. A technique as simple 
as the one given here can be improved in various ways; in doing this it 
should always be remembered that any improvement obtained at the cost 
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of much increased computation time will reduce the usefulness of the 
whole automatic reconstruction idea. The present edge detection scheme 
could readily be adopted for other classes of objects, e.g., those 
characterized by a sharp phase change or by a difference in texture. 

Once different values have been assigned to different types of texture 
it is possible to operate with the same edge detection routine on the 
thus interpreted array. 

In this investigation the spatial frequency approach to diffraction 
has been preferred to programming the Fresnel equation, mainly because 
repeated small propagation steps are possible when the transfer function 
is used, but also because it circumvents the Fresnel approximations and 
therefore gives better results in its region of validity. It was found 
that the necessary neglect of evanescent waves can give rise to very 
strong deterioration of the output, most of all for diffusely emitting' 
objects. It is tempting to use a lower sampling rate than the Nyquist 
rate because this avoids the necessity of neglecting evanescent waves: 
perfect reconstructions can be obtained this way, but only for computer 
generated diffraction patterns, or those obtained in an actual data 
acquisition system from periodic objects. For all other diffraction 
patterns the aliasing introduced by undersampling could be expected to 
be quite severe. 

For cases where the neglect of evanescent waves destroys edge 
information to a high degree, the reconstruction procedure developed 
in this investigation will not succeed since the resulting diffraction 
pattern at the object location does not belong to that class of 
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diffraction patterns which the computer accepts as objects because of 
the deterioration of the edge features. However, the test results 
presented are quite encouraging as they show that for many situations 
the computer can indeed make correct decisions as to the object location. 
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APPENDIX 



Main Loop of Subroutine ADERIV 

This computes the edgemeasure of the array AMAG 
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LISTING OF COMPUTER PROGRAM 
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SUBROUTINE SHUFL REARRANGES THE 4 QUADRANTS OF THE ARRAY DATA 
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